#include "datasourcegdal.h"
#include "mainwindow.h"

#include <gdal_priv.h>
#include <QDebug>
#include <osgDB/WriteFile>

#include <qmath.h>

#define IMGSIZE 256

rxtiler::DataSourceGDAL::DataSourceGDAL(const QString filePath)
{
    _xResolution=0.0;
    _yResolution=0.0;
    GDALAllRegister();//注册驱动
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");//支持中文路径
    poDataset = ( GDALDataset* )GDALOpen( filePath.toStdString().c_str(), GA_ReadOnly );

    double adfGeoTransform[6];
    CPLErr res = poDataset->GetGeoTransform(adfGeoTransform);
    if (res == CE_None)
    {
        _xResolution=adfGeoTransform[1];
        _yResolution=-adfGeoTransform[5];
    }
    dem = new osg::HeightField();

    int src_bandsize = poDataset->GetRasterCount();
    if(src_bandsize < 1){
                GDALClose(poDataset);
                poDataset = 0;
                return;
    }

    static int bandmap[3] = {1,1,1};
    static int bandmap3[3] = {1,2,3};
    static double nodata[3] = {0,0,0};

    if(src_bandsize < 3){
        GDALRasterBand* pband = poDataset->GetRasterBand(1);
        double nd = pband->GetNoDataValue();
        nodata[0] = nd;
        nodata[1] = nd;
        nodata[2] = nd;
        this->_bandmap = bandmap;
    }
    else{
        for(int i = 0; i < 3; i++){
            GDALRasterBand* pband = poDataset->GetRasterBand(i+1);
            double nd = pband->GetNoDataValue();
            nodata[i] = nd;
        }
        this->_bandmap = bandmap3;
    }
    this->_nodata = nodata;
}

rxtiler::DataSourceGDAL::~DataSourceGDAL()
{
    delete poDataset;
    poDataset = NULL;

}



rxtiler::Extent rxtiler::DataSourceGDAL::getExtent()
{
    double adfGeoTransform[6];
    if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None){
        int x=poDataset->GetRasterXSize();
        int y=poDataset->GetRasterYSize();
        double left=adfGeoTransform[0];
        double top=adfGeoTransform[3];
        double right=adfGeoTransform[0]+x*adfGeoTransform[1]+y*adfGeoTransform[2];
        double bottom=adfGeoTransform[3]+x*adfGeoTransform[4]+y*adfGeoTransform[5];
        _foemerTif.left = left;
        _foemerTif.right = right;
        _foemerTif.bottom = bottom;
        _foemerTif.top = top;


        return rxtiler::Extent(left,right,bottom,top);
    }else
    {
        return rxtiler::Extent(0,0,0,0);
    }
}

bool rxtiler::DataSourceGDAL::computerLevels(rxtiler::TileOption *option, int *min, int* max)
{
    rxtiler::Extent extent=getExtent();
    double deltX=extent.right-extent.left;//180
    double deltY=extent.top-extent.bottom;//360
    double delt=(deltX<deltY)?deltX:deltY;//180

    double resolution=(_xResolution<_yResolution)?_xResolution:_yResolution;

    if(delt>=36)
    {
        *min=0;
    }
    else
    {
        for(int i=0;i<21;i++)
        {
           double t= 36.0/(pow(2,i));
           if(t>=delt)
           {
               continue;
           }
           else
           {
               *min=i-1;
               break;
           }
        }
    }
    for(int i=0;i<21;i++)
    {
       double t= double(36.0/pow(2,i))/IMGSIZE;
       if(t>=resolution)
       {
           continue;
       }
       else
       {
            *max=i-1;
            break;
       }
    }

    return true;
}

osg::Image *rxtiler::DataSourceGDAL::getImage(rxtiler::TileKey &key, rxtiler::Extent &extent)
{
    if(!poDataset || !_tif.cover(extent)){
        return 0;
    }
    int out_bandsize = 3;
    int *out_bandmap = this->_bandmap;
    double* nodata = this->_nodata;

    int iWidth = poDataset->GetRasterXSize();
    int iHeight = poDataset->GetRasterYSize();

    //包含关系
    if (extent.compara(_tif) == -1)
    {
        int x = iWidth * (extent.left - _tif.left) / (_tif.right - _tif.left);
        int y = iHeight * ( _tif.top - extent.top) /  (_tif.top - _tif.bottom);
        int XSize = iWidth * extent.width() /  _tif.width();
        int YSize = iHeight * extent.height()  /  _tif.height();


        unsigned char *buffer = (unsigned char*)malloc(IMGSIZE*IMGSIZE*out_bandsize);
        CPLErr err = poDataset->RasterIO(GF_Read, x, y, XSize, YSize, buffer, IMGSIZE
                            , IMGSIZE, GDT_Byte, out_bandsize,out_bandmap,out_bandsize,IMGSIZE*out_bandsize,1);
        osg::Image* img = new osg::Image();
        img->allocateImage(IMGSIZE,IMGSIZE,out_bandsize,GL_RGBA,GL_UNSIGNED_BYTE);
        for(int x =0; x< IMGSIZE;x++){
            for(int y = 0; y< IMGSIZE; y++){
                double r = buffer[x*out_bandsize+y*IMGSIZE*out_bandsize] / 255.0;
                double g = buffer[x*out_bandsize+y*IMGSIZE*out_bandsize+1] / 255.0;
                double b = buffer[x*out_bandsize+y*IMGSIZE*out_bandsize+2] / 255.0;
                double a = 1;
                if(r == nodata[0] && g == nodata[1] && b == nodata[2] ){
                    r = g = b = a = 0;
                }
                img->setColor(osg::Vec4(r,g,b,a),x,IMGSIZE-1-y);
            }
        }
        free(buffer);

        return img;
   }

    //有重叠的关系
    else
    {
        //计算重叠区域
        overLapTif.right = _tif.right > extent.right ? extent.right : _tif.right;
        overLapTif.left = _tif.left > extent.left ? _tif.left : extent.left;
        overLapTif.bottom = _tif.bottom > extent.bottom ? _tif.bottom : extent.bottom;
        overLapTif.top = _tif.top > extent.top ? extent.top : _tif.top;
        //计算重叠区域在tif的位置
        int x = iWidth * (overLapTif.left - _tif.left) / (_tif.right - _tif.left);
        int y = iHeight * ( _tif.top - overLapTif.top) /  (_tif.top - _tif.bottom);
        int XSize = iWidth * overLapTif.width()/_tif.width()+1;
        int YSize = iHeight * overLapTif.height()/_tif.height()+1;
        //计算重叠区域在切片的位置
        int qX = IMGSIZE * (overLapTif.left - extent.left)/( extent.right - extent.left );
        int qY = IMGSIZE * ( extent.top - overLapTif.top)/(extent.top - extent.bottom);
        int qXSize = IMGSIZE * overLapTif.width()/extent.width()+1;
        int qYSize = IMGSIZE * overLapTif.height()/extent.height()+1;
        //切切切
        char *buffer = (char*)malloc(qXSize*qYSize*out_bandsize);
        poDataset->RasterIO(GF_Read, x, y, XSize, YSize, buffer, qXSize
                            , qYSize, GDT_Byte, out_bandsize,out_bandmap,out_bandsize,qXSize*out_bandsize,1);
        osg::Image* img = new osg::Image();
        img->allocateImage(IMGSIZE,IMGSIZE,4,GL_RGBA,GL_UNSIGNED_BYTE);
        int xmax = qX+qXSize < IMGSIZE ? qX+qXSize : IMGSIZE;
        int ymax = qY+qYSize < IMGSIZE ? qY+qYSize : IMGSIZE;
        for(int x =qX; x< xmax;x++)
        {
            for(int y = qY; y< ymax; y++)
            {
                int zx = x-qX;
                int zy = y-qY;
                double r = buffer[zx*out_bandsize+zy*qXSize*out_bandsize] / 255.0;
                double g = buffer[zx*out_bandsize+zy*qXSize*out_bandsize+1] / 255.0;
                double b = buffer[zx*out_bandsize+zy*qXSize*out_bandsize+2] / 255.0;
                double a = 1;
                if(r == nodata[0] && g == nodata[1] && b == nodata[2] ){
                    r = g = b = a = 0;
                }
                img->setColor(osg::Vec4(r,g,b,a),x,IMGSIZE-1-y);
            }
        }
        free(buffer);
        return img;
    }
}

osg::HeightField *rxtiler::DataSourceGDAL::getDEM(rxtiler::TileKey &key, rxtiler::Extent &extent)
{
    static int bandmap[1] = {1};;
    int iWidth = poDataset->GetRasterXSize();
    int iHeight = poDataset->GetRasterYSize();

    if (extent.compara(_tif) == -1)
    {
        int x = iWidth * (extent.left - _tif.left) / (_tif.right - _tif.left);
        int y = iHeight * ( _tif.top - extent.top) /  (_tif.top - _tif.bottom);
        int XSize = iWidth * extent.width() /  _tif.width();
        int YSize = iHeight * extent.height()  /  _tif.height();


        unsigned char *buffer = (unsigned char*)malloc(60*60);
        poDataset->RasterIO(GF_Read, x, y, XSize, YSize, buffer, 60
                            , 60, GDT_Byte, 1,0,0,0,0);

        dem->allocate(60,60);
        for(int x =0; x< 60;x++) {
            for(int y = 0; y< 60; y++){
                char h;
                h=buffer[60*y+x];
                dem->setHeight(x,y,h);
            }
        }

        free(buffer);
        return dem;
   }
   else
   {
        //计算重叠区域
        overLapTif.right = _tif.right > extent.right ? extent.right : _tif.right;
        overLapTif.left = _tif.left > extent.left ? _tif.left : extent.left;
        overLapTif.bottom = _tif.bottom > extent.bottom ? _tif.bottom : extent.bottom;
        overLapTif.top = _tif.top > extent.top ? extent.top : _tif.top;
        //计算重叠区域在tif的位置
        int x = iWidth * (overLapTif.left - _tif.left) / (_tif.right - _tif.left);
        int y = iHeight * ( _tif.top - overLapTif.top) /  (_tif.top - _tif.bottom);
        int XSize = iWidth * overLapTif.width()/_tif.width();
        int YSize = iHeight * overLapTif.height()/_tif.height();
        //计算重叠区域在切片的位置
        int qX = 60 * (overLapTif.left - extent.left)/( extent.right - extent.left );
        int qY = 60 * ( extent.top - overLapTif.top)/(extent.top - extent.bottom);
        int qXSize = 60 * overLapTif.width()/extent.width();
        int qYSize = 60 * overLapTif.height()/extent.height();
        //切切切
        unsigned char *buffer = (unsigned char*)malloc(60*60);
        poDataset->RasterIO(GF_Read, x, y, XSize, YSize, buffer, qXSize
                            , qYSize, GDT_Byte, 1,bandmap,1,qXSize,1);

        dem->allocate(60,60);
        for(int x =0; x< 60;x++) {
            for(int y = 0; y< 60; y++){
                if( (x >=qX && x <= qX+qXSize) && ( y >= qY && y<= qY+qYSize ) )
                {
                    int zx = x-qX;
                    int zy = y-qY;
                    float h = buffer[zx+zy*qXSize];
                    dem->setHeight(x,y,h);
                }
                else
                {
                    dem->setHeight(x,y,0);
                }
            }
        }

        free(buffer);
        return dem;
    }

}
#include "GeoRaster.h"

osg::HeightField *rxtiler::DataSourceGDAL::getDEM2(rxtiler::TileKey &key, rxtiler::Extent &extent)
{
    int dem_size = 60;
    if(!_tif.cover(extent)){
        osg::ref_ptr<osg::HeightField> dem1 = new osg::HeightField();
        dem1->allocate(dem_size,dem_size);
        return dem1.release();
    }

    int iWidth = poDataset->GetRasterXSize();
    int iHeight = poDataset->GetRasterYSize();


    Extent ext2 = extent.intersect(_tif);
    int x = iWidth * (ext2.left - _tif.left) / _tif.width();
    int y = iHeight * ( _tif.top - ext2.top) / _tif.height();
    int XSize = ::ceil(iWidth * ext2.width() /  _tif.width());
    int YSize = ::ceil(iHeight * ext2.height()  /  _tif.height());

    ext2.left = _tif.left + _tif.width() * x / iWidth;
    ext2.right = _tif.left + _tif.width() * (x+XSize) / iWidth;

    ext2.top = _tif.top - _tif.height() * y / iHeight;
    ext2.bottom = _tif.top- _tif.height() * (y+YSize) / iHeight;
    int bxSize = ::ceil(ext2.width() * dem_size / extent.width());
    int bySize = ::ceil(ext2.height() * dem_size /extent.height());
    osg::ref_ptr<osg::HeightField> dem1 = new osg::HeightField();
    dem1->allocate(bxSize,bySize);
    if(bxSize <=1 || bySize <=1){
        return dem1.release();
    }
    int xiuzhengx = XSize / bxSize + 0.5;
    int xiuzhengy = YSize / bySize + 0.5;
    int l = x - xiuzhengx * 0.5;
    int t = y - xiuzhengy * 0.5;
    int r = x + XSize + xiuzhengx;
    int b = y + YSize + xiuzhengy;
    if(l < 0) l = 0;
    if(t < 0) t = 0;
    if(r > iWidth) r = iWidth;
    if(b > iHeight) b = iHeight;


    short* buffer = (short*)malloc(bxSize * bySize * 2);
    poDataset->RasterIO(GF_Read,l,t,r-l,b-t,buffer,bxSize,bySize,GDT_Int16,1,0,0,0,0);

    for(int x =0; x< bxSize;x++) {
        for(int y = 0; y< bySize; y++){
            short h;
            h=buffer[bxSize*y+x];
            dem1->setHeight(x,y,h);
        }
    }
    free(buffer);
    if(ext2.equal(extent) && bxSize == dem_size && bySize == dem_size){
        return dem1.release();
    }
    osg::HeightField* dem2 = new osg::HeightField();
    dem2->allocate(60,60);
    mapRaster(dem1,ext2,dem2,extent);

    return dem2;
}

#include "geoimage.h"

osg::Image *DataSourceGDAL::getImage2(TileKey &key, Extent &extent)
{
    if(!poDataset){
        return 0;
    }
    int img_size = IMGSIZE;
    if(!extent.cover(_tif)){
        osg::Image* img = new osg::Image();
        img->allocateImage(img_size,img_size,4,GL_RGBA,GL_UNSIGNED_BYTE);
        return img;
    }
    if(_tif.contain(extent)){
        return _getContianImage(key,extent);
    }
    int out_bandsize = 3;
    int *out_bandmap = this->_bandmap;
    double* nodata = this->_nodata;

    int iWidth = poDataset->GetRasterXSize();
    int iHeight = poDataset->GetRasterYSize();


    Extent ext2 = extent.intersect(_tif);
    int x = iWidth * (ext2.left - _tif.left) / _tif.width();
    int y = iHeight * ( _tif.top - ext2.top) / _tif.height();
    int XSize = ::ceil(iWidth * ext2.width() /  _tif.width());
    int YSize = ::ceil(iHeight * ext2.height()  /  _tif.height());

    int bxSize = ::ceil(ext2.width() * img_size / extent.width());
    int bySize = ::ceil(ext2.height() * img_size /extent.height());
    if(bxSize <=1 || bySize <=1){
        osg::Image* img = new osg::Image();
        img->allocateImage(img_size,img_size,4,GL_RGBA,GL_UNSIGNED_BYTE);
        return img;
    }

    char* buffer = (char*)malloc(bxSize * bySize * out_bandsize);
    poDataset->RasterIO(GF_Read, x, y, XSize, YSize, buffer, bxSize
                        , bySize, GDT_Byte, out_bandsize,out_bandmap,out_bandsize,bxSize*out_bandsize,1);
    osg::Image* img = new osg::Image();
    img->allocateImage(bxSize,bySize,4,GL_RGBA,GL_UNSIGNED_BYTE);
    for(int x =0; x< bxSize;x++) {
        for(int y = 0; y< bySize; y++){
            int idx = (bxSize * y + x)*out_bandsize;
            double r = buffer[idx] / 255.0;
            double g = buffer[idx+1] / 255.0;
            double b = buffer[idx+2] / 255.0;
            double a = 1;
            if(r == nodata[0] && g == nodata[1] && b == nodata[2] ){
                r = g = b = a = 0;
            }
            img->setColor(osg::Vec4(r,g,b,a),x,bySize-1-y);
        }
    }
    free(buffer);

    if(ext2.equal(extent) && bxSize == img_size && bySize == img_size){
        return img;
    }
    osg::Image* img2 = new osg::Image();
    img2->allocateImage(img_size,img_size,4,GL_RGBA,GL_UNSIGNED_BYTE);
    mapImage(img,ext2,img2,extent);
    return img2;
}

void rxtiler::DataSourceGDAL::setNowLrbt(const rxtiler::Extent &lrbt)
{
    _tif.left = lrbt.left;
    _tif.right = lrbt.right;
    _tif.top = lrbt.top;
    _tif.bottom = lrbt.bottom;
}

rxtiler::Extent rxtiler::DataSourceGDAL::getFoemerLrbt()
{
    return _foemerTif;
}

osg::Image *DataSourceGDAL::_getContianImage(TileKey &key, Extent &extent){
    int out_bandsize = 3;
    int *out_bandmap = this->_bandmap;
    double* nodata = this->_nodata;

    int iWidth = poDataset->GetRasterXSize();
    int iHeight = poDataset->GetRasterYSize();
    int x = iWidth * (extent.left - _tif.left) / (_tif.right - _tif.left);
    int y = iHeight * ( _tif.top - extent.top) /  (_tif.top - _tif.bottom);
    int XSize = iWidth * extent.width() /  _tif.width();
    int YSize = iHeight * extent.height()  /  _tif.height();


    unsigned char *buffer = (unsigned char*)malloc(IMGSIZE*IMGSIZE*out_bandsize);
    CPLErr err = poDataset->RasterIO(GF_Read, x, y, XSize, YSize, buffer, IMGSIZE
                        , IMGSIZE, GDT_Byte, out_bandsize,out_bandmap,out_bandsize,IMGSIZE*out_bandsize,1);
    osg::Image* img = new osg::Image();
    img->allocateImage(IMGSIZE,IMGSIZE,out_bandsize,GL_RGBA,GL_UNSIGNED_BYTE);
    for(int x =0; x< IMGSIZE;x++){
        for(int y = 0; y< IMGSIZE; y++){
            double r = buffer[x*out_bandsize+y*IMGSIZE*out_bandsize] / 255.0;
            double g = buffer[x*out_bandsize+y*IMGSIZE*out_bandsize+1] / 255.0;
            double b = buffer[x*out_bandsize+y*IMGSIZE*out_bandsize+2] / 255.0;
            double a = 1;
            if(r == nodata[0] && g == nodata[1] && b == nodata[2] ){
                r = g = b = a = 0;
            }
            img->setColor(osg::Vec4(r,g,b,a),x,IMGSIZE-1-y);
        }
    }
    free(buffer);

    return img;
}
